Two interacting fermions in a ID harmonic trap: matching the Bethe ansatz and 

variational approaches 
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In this work, combining the Bethe ansatz approach with the variational principle, we calculate 
the ground state energy of the relative motion of a system of two fermions with spin up and down 
interacting via a delta-function potential in a ID harmonic trap. Our results show good agreement 
with the analytical solution of the problem, and provide a starting point for the investigation of 
more complex few-body systems where no exact theoretical solution is available. 



1. INTRODUCTION 

Few-body quantum systems composed of atoms and 
molecules are some of the most simple structures that 
constitute the building blocks of matter. Despite this 
simplicity their study has recurrently been challenging. 
One of the reasons being that one cannot make use of 
standard statistical methods and still there are enough 
degrees of freedom to make it a complex problem, often 
not solvable for as few as three bodies. The interest 
in few body systems is manifold and has over the time 
appeared in nuclear and particle physics as well as in 
atomic and molecular studies. Currently, a renewed in- 
terest also emerged in relation to the experimental study 
of Bose and Fermi gases, since few body interactions 
may play a far from trivial role in their behavior 
Moreover,the recent and impressive development of the 
technology associated to the study of the Bose-Einstein 
condensation phenomena in fields such as ultracold 
gases, Mott insulators and optical lattices led to the 
possibility of controlling in an increasingly precise way 
the number of atoms trapped in a well. 

In particular, a great deal of interest has been devoted 
to the study of distinguishable trapped few-fermion sys- 
tems. The most recent experimental achievement being 
the realization of a system of two fermionic atoms of 
6 Li with tunable interactions 0, In this experiment, 
the ground state energy of the system was measured 
and compared to an analytical result that exists in this 
particular case [J] (see also Q), which however is not 
extendable if one includes more atoms. Therefore a good 
approximation that may be generalized to more than 
two atoms is of interest. It is worthwhile to mention 
here that the Hamiltonian employed in the calculation 
of the ground-state energy of this 2-fermion experiment 
is basically equivalent to the one used to discuss the 
existence of exotic pairing mechanisms closely related 
to the elusive Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) 
[H, 0] state. In that case one deals with a higher number 
of particles, in addition to a spin imbalance, and the 
thermodynamical Bethe ansatz coupled to a local density 
approximation was used to discuss the resulting phase 
diagrams and density profiles of trapped fermionic 6 Li 
atoms in ID tubes [a-UOj. 



With the above motivations, viewing the prospect 
of new few body experiments [ll|, we propose an 
alternative possibility, a variational approach based on 
the use of the Bethe-ansatz solution for a system with 
delta-function interactions. Our choice will take into 
account the knowledge of the exact solution of the one 
dimensional many-body system with repulsive or attrac- 
tive delta-function potentials fl2l - fl4) and consider the 
trapping as a kind of perturbation. By this we mean that 
the bulk of our ansatz is supposed to grasp the behavior 
of the interacting particles which happen to be trapped 
in a harmonic well. To the extent of our knowledge this 
is an unexplored possibility and, for that matter, one 
that has the potential to be systematically generalized 
from two to more particles. Our approach consists in 
calculating the ground state of the few-fermion model 
having in mind the variational principle, such that the 
actual ground state energy is smaller than the ground 
state of the Hamiltonian with delta interactions, which 
we know exactly by the Bethe-ansatz methods, plus a 
part that is the mean value of the harmonic potential 
for our ansatz. 

In the following we develop our systematics for the 
variational calculation of the ground state of a two- 
fermion system. The next section will be devoted to set 
forth the system of two interacting fermions that we are 
investigating, then in section 3 we introduce our varia- 
tional ansatz, which as mentioned is inspired in a paradig- 
matic solution for one dimensional systems |12h15| . the 
Bethe-ansatz [l^. In section 4 we present our results for 
the repulsive and attractive cases and in section 5 these 
results are compared with the one obtained in 0, HJ for 
the relative motion. In the Appendix we provide details 
concerning the construction of the Bethe ansatz part in 
absolute coordinates and briefly discuss its extension to 
the general case of N fermions. 



2. SYSTEM 

Let us consider a system of two interacting fermions, 
for instance two fermionic atoms with mass m, in an axi- 
ally symmetric harmonic trap with angular frequency u>. 
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Such a system can be described by the following Hamil- 
tonian: 



h 2 d 2 h 2 d 2 
H = -^dx{-^dxl + VA{x ^ X2) + Vl{Xl - X2h 

where x\ and x 2 denote the position of the two fermions 
and Va (xi, x%) represents the trapping potential. 



ansatz concept dMHiil. The parameter a controls the 
decay of the trial function outside the trap and L indi- 
cates the limit where this decay starts. Inside the trap, 
where the contact interaction is relevant the trial func- 
tion will take the form of the Bethe ansatz. As usual a 
variation of these parameters provides a minimal value, 
which should approximate the ground state of the sys- 
tem. In that way we have a wavefunction that gives a 
realistic picture of the physical processes involved. 



Va {xi, x 2 ) — ^moj 2 xi 2 + ^mw 2 i 2 2 - (2.2) 

For sufficiently low energies the interaction potential Vj 
can be taken as a delta-function contact potential, such 
that, 



Vi (xi - x 2 ) = 2cS (x 2 - xi) , 



(2.3) 



where c is the interaction strength. The potential is 
repulsive or attractive, respectively for c > or c < 0. 

Here, as the harmonic potential and the kinetic energy 
are quadratic it is convenient to separate the relative mo- 
tion from the center of mass motion. This can be easily 
attained by using center of mass and relative coordinates 
given by, 



v xi+x 2 
Jt = ~ , x 



x 2 - XX- 



(2.4) 



One can then decompose the total Hamiltonian in the 
center of mass Hqm and relative motion H re i parts, 



3. ANSATZ 

In the present section we exhibit our variational 
ansatz. Further details concerning the construction of 
the Bethe ansatz part can be found in the Appendix, 
where we also briefly discuss its extension. 

As shown in Fig. 1 there are three relevant regions for 
a wavefunction of our system. We can delineate these 
regions by the parameter L. Our variational ansatz as- 
sumes then the following configuration: 



lp! = e -«(*+ L ) i/jjj (-L) 



—00 < x < —L 



i>u = (e 



ikL „—ikx 



^ikx\r\ 



ikL „ikx 



-ikx\Q 



)Q(x) 

)6(-x) -L < x < +L 



where O is the Heaviside step function. 



+L < x < +00 
(3.1) 



H, 



d 2 



1 



CM 



2M dX 2 



-Muj z X 



1 v 2 



(2.5) 



|^+2c 5 (,) + ^V. (2.6) 



In the above, /1 = ?r is the reduced mass and M = 2m 
the total mass. It can be seen that the eigenfunctions 
and eigenenergies of Hcm are those of the harmonic 
oscillator. Notice that for the general N case, by the use 
of Jacobi coordinates, the Hamiltonian is also separable 
(see Appendix). 




Region I 



Region II 



Region III 



Now, for the Hamiltonian H re i , we shall apply the vari- 
ational principle, 



Egs < 



i\ H re i \ tp) 



(2.7) 



where tp (x) is a continuous trial wavefunction. 



Figure 1: Schematic representation of the normalized proba- 
bility density, ^(rr)! 2 in the relative coordinates system. The 
variational parameter L may be used to delineate three re- 
gions according to the boundaries with respect to the har- 
monic potential. 



The novelty in our approach is that the trial function, 
which we denote ip (x, {a, L}), will encompass the Bethe 



In region II {—L < x < +L) the two fermion system 
is subject to the contact potential and the harmonic 
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trap. Due to the symmetry of the system, in the vicinity 
of the central axis, where (x — 0), the interaction term 
is dominant. In other words, any contact interaction 
takes over the harmonic potential. For this reason we 
assume that is possible to approximate the wavefunc- 
tion in region II by the wavefunction that describes 
a system with two distinct fermions with a contact 
interaction. Historically, such systems where studied 
in one dimensional lattices of size "L" and periodical 
boundary conditions, being exactly solved in [l3|, Il4j . 
Later major contributions for this problem where given 
by [Hi followed by others, such as [l7| and [l8l-|2ll|. 

Therefore, our choice for a trial wavefunction in this 
region corresponds to the Bethe ansatz solution for 
fermions interacting trough a delta function potential in 
relative coordinates, such that ipxi is built as the eigen- 
function of the interaction Hamiltonian 



H int 1p 



II 



n 2 d 2 

W + M(a:) . 



ipn = —k tpn. 
2\i 



More importantly, this means that ipu will correspond 
to the eigenfunction of Hint with energy E{ n t = f^fc 2 for 
all the "quasi-momenta" k's that satisfy the following 
equation 



ikL 



k + i f C t 

k — ijjpc 



(3.2) 



known as the Bethe ansatz equation. 



Before we proceed to the analysis of the repulsive and 
attractive cases we have still to explain how to deal with 
the continuity of the wavefunction on the boundaries be- 
tween the three regions. The continuity condition in all 
the interval dictates that 

"*Pi (-L) = iPii(-L), 
il>u{+L) = tp UI (+L). 

In regions / and III the harmonic potential is the only 
one present, so the simplest choice that take into account 
the system behavior should be an eigenfunction of the 
harmonic oscillator Hamiltonian. As we expect a rapid 
decay of the probability density in these regions, ipi and 
tpin have the form of a Gaussian and depend on an- 
other variational parameter, a. It is important to notice 
that although our choice is continuous for all x, its first 
derivative is not. Later we also consider the contribution 
of this discontinuity to the total value of the ground state 
energy. 



4. RESULTS 

4.1. Repulsive interaction, (c > 0) 

In this case just some values of k, purely real ones, 
satisfy (j3.2p for the ground state, therefore in this sub- 
section we only consider k € K. In order to apply the 
variational principle, we first compute the normalization 
of the wavefunction, which yields 



Hence, for each value of the coupling c we need 
to determine which value of the quasi-momentum k 
satisfies (|3.2p . for the ground state, in order to have ipn 
completely defined. 



A careful analysis of p.2[) shows that the possible 
values for the quasi-momenta k depend on the sign of 
c fl2| . For the repulsive case (c > 0), only real fc's are 
ground-state solutions of the Bethe ansatz equations 
(13. 2p and, accordingly, are the values entering in ipn. 
For the attractive case (c < 0) the fc's composing the 
ground state are pure imaginary numbers. We will then 
consider both cases separately. 



ip) = 2J^- [1 + cos (kL)] + y [kL + sin (kL)] (4.1) 
V 2a k 



and then the expectation value of H re i , which value is 



h 2 a flT rv , 2h 2 k 2 L 
H V 2a /i 



+ /Ltcj [l + cos (kL)] 




>;3 [2L 3 fc 3 + 3(fc 2 L 2 -l)sin(fcL) + 3fcLcos(fcL)] . 

(4.2) 



Our proposal consists then in building the central part 
of our trial function as the Bethe ansatz wavefunction for 
the relative motion of two distinct fermions interacting 
via a delta function. This problem is completely solvable 
and in its generality applied to any number of fermions 
[l2TJl4l |. Notice that in the literature, one usually 
considers the contact interaction as a perturbation to 
the harmonic potential Hamiltonian. We show here 
how to use the full strength of the Bethe ansatz in a 
variational approach. 



In the expression (|4.2p we also considered the contri- 
bution from the discontinuity of the wavefunction first 
derivative at x = ±L, that is: 



lim 



-L+t rL+e 

1p* (H re itjj) dx+ tfj* (Hreit/j) dx 



—L— 

2h 2 k 



L-£ 



sin (kL) , 



which comes from the kinetic term of H 



rcl • 
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The Bethe ansatz equations (I3.2[) for the ground state 
in the repulsive interaction reduce then to 



In order to apply the variational principle, we first com- 
pute the normalization of the wavefunction, which yields 



2 (2[ic 
L aiCtan {Wk 



which are much simpler to solve. 



(4.3) 



We have then all the necessary ingredients to proceed 
with the numerical minimization of ^^p^]^ with 
respect to the parameters a and L. Basically to each 
assigned c, we sweep over all values of a and L, calculate 
k for each L and establish the parameters a* and L* 
such that ^^j^j^ takes the least possible value. In 
this way we determine the ground state energy of the 
two fermion system as a function of the coupling c via 
the variational principle, where the trial wavefunction is 
constructed by means of the Bethe ansatz. This result 
is depicted in Fig. 2 using the physical variables e and 
o>id- We give more details in the next section where we 
also compare this result with the analytical solution (ij. 



Limiting case: Notice that in the limit L — > 0, c - 
(harmonic oscillator) the expression (|4.2p reduces to 







(V| V) = 2e 



-k'L 



' — [1 + cosh (k'L)] 
2a 



- [Lk' + sinh (k'L)} 



and then the mean value of H re i , obtaining 
(VI H rel |V) = —\j]^e- k ' L [1 + cosh(fc'L)] 



2h 2 k' 2 L^_ k , L + ^_ k , L 3fc/: 



(4.8) 



-3 (k' 2 L 2 + l) sinh (k'L) - 3k' L cosh (k'L 
+ Hu 2 e- k ' L [1 + cosh (k'L)} ' 



where again we considered the contribution of the dis- 
continuity of the wavefunction derivative at the points 
x = ±L, 




(4>\H rel |V) h 2 a fiuj 2 

iim ; ; : = 1 

l^o (V| V) 2^ 8a 



(4.4) 



Upon extremization of the total energy with respect to 
a in the limit L — > 0, the minimum value is the one for 
the value a* of the parameter 



lim 



V* (H re ltp) dx ■ 



2h 2 k ' 



-k'L 



sinh (kL) , 



into the expression (|4.9|) . 



L+e 



V* (Hrellp) dx 



L-t 



such that 



Hrel IV) 



lim 

L^O (V| V) 



= 2^' 



(4.5) 



(4.6) 



which, as expected, is simply the ground state energy of 
the harmonic oscillator. 



4.2. Attractive interaction, (c < 0) 

In this case only the purely imaginary values of k sat- 
isfy the Bethe ansatz equations (|3.2p for the ground state; 
for this reason we will consider k G C. Thus, it is con- 



venient to define k — ik! , kl £ 
write (13.21) as 



In terms of k' we can 



-k'L 



k' 



2M, 



(4.7) 



which can be solved by numerical methods. 



As in the attractive case we numerically minimize the 
(ip\H rcl \ip) w ^ n reS p ec |; t th e parameters 



mean energy 



WW) 



a and L and determine the ground state energy of the 
system as a function of the coupling c. This result is 
depicted in Fig. 2. In the next section we make a com- 
parison with the analytical solution obtained in Q ■ 

Limiting case: Notice again that in the limit L —> 
0, c —> (harmonic oscillator) 



lim ■ 



(V| V) 



io a uuj 

^ + ^- 4.10 
2/x 8a 



Upon extremization of the total energy with respect 
to a in the limit L — >• 0, we obtain |a*| = such that 
we find again the ground state energy of the harmonic 
oscillator. 



5. COMPARISON 

The results obtained in the previous section by means 
of the Bethe ansatz and the variational principle to find 
the two fermion system ground state as a function of the 
coupling parameter are presented in Fig. 2 (red straight 
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line) . For convenience and to compare with known re- 
sults we are using the following variables, 



Egs 

fkd 



aiD 



ft 3 / 2 
2c Y M 



(5.1) 



where e denotes the energy of the ground state shifted 
by the zero point energy in hw units and am is the one 
dimensional scattering length. 



The analytical solution in relative coordinates for a 
system of two distinct fermions interacting via a delta 
function and confined in a harmonic trap was first 
obtained in Basically in this work they expanded 
the unknown wavefunction in a complete set of the 
simple harmonic oscillator functions. Later, these results 
were generalized to different geometries of the trapping 
potential in [f|, and among others in [22|-[26[. 



a 1d 



Figure 2: Energies for the ground state of the relative mo- 
tion for two distinct fermions interacting via a delta function 
potential and confined in a harmonic trap of frequency w. 
The analytically obtained energy levels (black dotted line) 
are compared to the results for the energy given by the com- 
bination of the Bethe ansatz with the variational principle 
(red straight line). 



CONCLUSION 



Essentially, the following implicit equation determin- 
ing the eigenenergies of the relative motion in a one- 
dimensional harmonic potential was obtained 



2o,D = fTTfF <5 ' 2) 



where T (x) is the complete gamma function. 



This solution for the ground state is plotted in Fig. 2 
(black dotted line). We can observe a very good agree- 
ment between this result and the result that we obtained 
combining the Bethe ansatz and the variational principle 
(red line). This places our ansatz as a potential candi- 
date for the extension to more than two fermions, were 
an analytical solution does not exist. The fact that the 
measured properties of this system Q may, with a good 
agreement, be compared with the theoretical results [4| 
makes this subject even more captivating. 



In this work we obtained the ground-state energy of 
two distinct fermions in a ID harmonic trap within a vari- 
ational approach, but from a distinct perspective, aim- 
ing a new view for the problem of few fermions. The 
reasoning beneath our variational ansatz choice was to 
exploit the exact solution for the one dimensional system 
of fermions interacting by means of a contact potential 
solution, the Bethe ansatz. Usually, in the literature, 
one takes a route different from ours by considering the 
harmonic trap Hamiltonian as the relevant one and the 
contact interaction as a perturbation. But, since for delta 
function interactions, we have at our disposal the Bethe 
ansatz technology it is almost natural to use it. Thus, 
we chose a trial wavefunction for this system that con- 
tains a great deal of information about the physics of the 
two fermions inside the trap and supplement it by the 
knowledge of the harmonic oscillator Hamiltonian. The 
good agreement between our results and existing analyt- 
ical results shows that our ansatz fulfills our expectation 
and has the potentiality to shed light on the spectrum 
of strongly correlated few-body quantum systems. Us- 
ing the methods established in this work it is in principle 
possible to extend our studies to more complex systems, 
composed of three or more fermions, which are currently 
of experimental interest [llj , and also in this case one can 
profit of the exact solution for the contact interaction. 
The procedure for higher N brings however a substan- 
tial operational growth as there are Nl x AH coefficients 
of the Bethe ansatz to be determined and the number 
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of regions of the complete variational ansatz, as in (|3.1I) . 
increases correspondingly. 
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how to apply the same principles for a higher number of 
fermions. 

The rationale we use in our construction is the Bethe 
ansatz method for obtaining the energy spectra of exactly 
solvable Hamiltonians. Let us then consider two fermions 
interacting through a delta function potential in a one 
dimensional system with periodicity L, which has the 
following Hamiltonian 



n 2 d 2 n 2 d 2 n .. 



Appendix 

We develop here, in detail, how we built the trial wave- 
function for the two body problem and then we indicate 



where xi and x 2 are the position of each fermion and 
c is the interaction strength. The most general wave- 
function for such system in absolute coordinates in the 
region X\, x 2 £ [-L/2, +L/2] can be written as [HJ 



ip(x 1 ,x 2 ) = \a{le tiklXl+k2X2) + alle^ k2Xl+klx ^] Q(x 2 - x x ) 

+ \al\e l{klX2+k2Xl) + a 2 2 \e l ^ k2X2+klXl A 6 {x x - x 2 ) ■ 



(A-2) 



where ki and k 2 are the "quasi-momenta" for the 
fermions and the coefficients "a" are to be determined 
by physical arguments. The action of the Hamiltonian 
on the wavefunction results in 



this implies that 



n 12 
"12 



'21 



21 

a 12 = a 12 

^21 — „ 
a 21 = a 21 



(A-6) 



Hip(xi,x 2 ) = Eip(xi,x 2 ) + undesirable terms, (A-3) Besides, the periodic boundary conditions 



where the "undesirable terms" are functions of k\ and 
k 2 . When, as usual, one requires that these terms be 
null, k\ and k 2 must satisfy certain consistency relations 
known as the Bethe ansatz equations [12]. These depend 
on the wavefunction symmetry Energy and momentum 
are given respectively by 



E 



k 2 > 



K = k-L + k 2 



(A-4) 



Once the system is in a spin singlet configuration (antis- 
symetric) the wavefunction must be spatially symmetric, 



ip (x 2l xi) = ip (xi,x 2 ) , 



(A-5) 



i>(x j =-L/2)=i>(x j = +L/2), j — 1,2 (A-7) 
lead to the relations 



ai2 = a 2 ie 
a,2i = a 12 e 



(A- 



It can be shown that, for the ground state, K = => 
k 2 = —k\. Therefore, we can write the symmetric wave- 
function with periodical boundary conditions in terms of 
the absolute coordinates as 
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ip(xi,x 2 ) = a 21 



ikiL^ikifxi— ra) _|_ g — ifei(a;i — X2) 



9 (£2 - Xi) 



(A-9) 



+a 21 Lifei^ e -^i(^i-^) + e «i(ai-« 3 )j e ( Xl _ X2 ) 



Then, for convenience, defining x = x\ — x 2 and fc = k\, 
we obtain our ansatz in the relative coordinates system: 

V> (x) = a 2 i [e lfc V fca: + e~ lkx ] 6 (-a;) (A-10) 
+a 21 [e ikL e~ ikx + e lkx ] S (x) 

The constant 021 is obtained by the normalization condi- 
tion. The above wavefunction is the eigenfunction of the 
interaction Hamiltonian in relative coordinates in the 
—L<x<L interval, and constitutes the central part of 
our ansatz. In the x — > ±00 limits it is expected that the 
wavefunction exhibits an asymptotic behavior similar to 
the behavior of an harmonic oscillator wavefunction with 
exponential decay. Then the trial wavefunction in the 
other intervals is obtained by the condition of continuity 



of the wavefunction on the boundaries between all 
regions. 



The generic N fermion case can be dealt with in a 
completely analogous way. The wavefunction in abso- 
lute coordinates is similar to that of Eq. (|A-2jl with co- 
efficients a ■■ (the indices run from 1 to N). It is 
possible then to proceed in exactly the same way as be- 
fore for all other terms. Requiring the same physical 
principles as above for all the N\ regions of the type 
%l < x 2 < ... < xn £ [0,L] respecting the Bethe ansatz 
with periodicity L, a complete wavefunction in absolute 
coordinates can be constructed. For instance, for the case 
N = 3 we have, in a compact form: 



3 3 

ip(x 1 ,x 2) x 3 )= J2 E ^:^ x ' +k ^ +k ^e( Xn -x m )e(x m -x l ) 

Z,m,n— 1 X,fi,i/— 1 



An important element here is that the coefficients of 
the Bethe ansatz, are related through a well established 
transformation where the operators are given by [T2I [27J 
and satisfy the Yang-Baxter equation (l3l . [l5j . It is im- 
portant to note that, in this case, it is convenient to 
change from the absolute coordinates system to Jacobian 
coordinates and it is possible to show that for both the 
contact and trapping interaction the resulting Hamilto- 
nian is separable in center of mass and relative coordi- 



nates (28J. This enable us to, knowing the Bethe ansatz 
result for absolute coordinates, obtain the result for the 
relative coordinates Hamiltonian. In other words the de- 
scribed procedure, when N is increased, though cumber- 
some (the number of coefficients increases as AH 2 ), allows 
one to build the Bethe ansatz part of the whole varia- 
tional ansatz. The caveat, of course, is that the number 
of regions, such as in (|3.1[) . where one has to use conti- 
nuity conditions also increases accordingly. 
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